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SUMMARY 

The  main  objective  of  this  study  is  to  develop  an  analytical  model 
for  the  prediction  of  the  mixing  and  decay  of  a  passive  scalar  within 
the  turbulent  wake  of  a  submarine  operating  beneath  the  ocean  surface. 

Even  though  the  interest  is  on  the  passive  scalar,  the  essence  of  the 
problem  is  the  Jetermi nation  of  the  dynamics  of  a  submarine  wake.  Once 
the  wake  flow  field  is  determined,  the  decay  of  the  concentration  of  a 
passive  scalar  can  then  be  easily  obtained. 

Depending  on  the  distance  benind  the  submarine,  two  somewhat 
different  technical  p«xblems  are  encountered  in  the  analytical  modeling 
of  the  wake.  For  the  region  relatively  close  to  the  submarine,  (wake  age 
of  the  order  of  an  hour),  the  flow  field  is  influenced  by  the  effects  directly 
caused  by  the  body.  The  Droper  analytical  model  must  take  into  consideration 
the  following  characteristics  of  a  submarine  wake:  the  turbulence 
generated  by  the  body  and  the  propeller;  the  momentum  less  wake  behind  a 
submarine  moving  at  a  constant  speed,  and  the  stratified  nature  of  the  fluid 
media.  For  the  region  far-  behind  the  submarine,  the  sub-generated 
perturbations  decay  below  the  natural  background  turbulence  levels. 

Further  diffusion  of  the  passive  scalar  in  this  region  of  the  wake  is 
expected  to  be  dominated  by  the  ocean  background  turbulence  and  the 
character  of  the  diffusion  plume. 

An  analytical  model  for  the  near  wake  region,  based  on  the  concepts 
cf  local  similarity  and  turbulent  entrainment  of  the  external  fluid,  is 
formulated  in  Section  2.  A  turbulent  model  is  introduced  which  pays 
particular  attention  to  the  effect  of  stratification.  Using  basically 
an  integral  conservation  approach,  the  problem  is  reduced  to  solving  a 
set  of  five  ordinary  differential  equations  with  given  initial  conditions. 

By  properly  normalizing  the  governing  equations,  two  primary  sca’ing 
parameters  for  the  submarine  wake  are  obtained.  The  most  important 
scaling  parameter  is  the  turbulent  Froude  number  Fr  =  u^/NL,  where  u^Q 
denotes  the  initial  turbulent  intensity,  N  =  (-S-^)1^2  the  Brunt-Vasaila 
frequency  and  L  is  the  initial  wake  width.  To  a  lesser  degree,  the  initial 
turbulent  intensity  level  u^J0/U0  provides  another  governing  parameter. 
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Theoretical  or  empirical  estimates  of  the  value  of  these  parameters,  for 
any  particular  experiment  or  field  condition  is  not  available  at  the 
present  time.  However,  a  parametric  study  will  soon  be  performed  and 
comparisons  with  existing  laboratory  experiments  will  be  attempted.  These 
results  will  be  discussed  in  the  final  report.  Meanwhile,  effort,  both 
analytically  and  experimentally,  should  be  directed  toward  a  logical  and 
consistent  determination  of  the  value  of  these  governing  parameters. 

The  effects  of  the  background  turbulent  ocean  dynamics  has  been 
completely  ignored  in  the  formulation  presented  in  Section  2,  which  is 
valid  fcr  the  relatively  near  wake  region.  As  a  first  step  toward 
increasing  the  understanding  of  this  important  effect,  an  In-depth  review 
is  presented  in  Section  3  of  existing  theoretical  and  experimental 
studies  of  diffusion  in  an  ocean  environment. 

In  analytically  modeling  the  diffusion  of  a  passive  scalar,  the 
diffusion  coefficients  K.  must  be  known  for  both  vertical  and  horizontal 
transport.  The  review  in  Section  3,  concludes  that  horizontal  diffusion 
in  the  ocean  dominates  over  vertical  diffusion  (Kh  >  104KV)  because  of 
the  effect  of  the  stratified  media  in  retarding  vertical  turbulence 
fluctuations.  Thus  diffusion  from  a  moving  continuous  source  proceeds 
as  a  thin  fan  in  the  opposite  direction  from  the  source  trajectory. 

The  horizontal  diffusion  coefficient  is  estimated  with  the  aid  of  the 
theories  of  Taylor,  Batchelor  and  Richardson,  showing  that  Kh  is  proportional 
to  the  4/3  power  of  the  scale  of  diffusion  i.  Examination  of  a  number 
of  dye  diffusion  experiments  in  the  ocean  indicates  that  is  approximately 
represented  by  the  Richardson  4/3  law,  (K  ^  £1,15)  and  that  the  lateral 
dimension  o  of  the  dye  plume  increases  rapidly  with  time  a  *  t1*17  (a  ^  t3/2 
for  the  Richardson  law).  This  rate  of  growth  indicates,  for  example,  that 
the  size  of  a  dye  patch  would  be  on  the  order  of  several  kilometers  after 
a  period  of  diffusion  of  one  day. 

In  Section  3,  the  diffusion  coefficients  and  rate  of  spread  are 
used  to  estimate  the  concentration  decay  behind  a  continuous  source. 

For  a  point  source,  it  is  shown  that  the  concentration  falls  off 
inversely  as  the  diffusion  time  to  the  1.17  power  (or  3/2  power  for  the 
Richardson  law).  Far  field  solutions  are  also  presented  for  diffusion 
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from  a  finite  bounded  concentration  field.  Such  initial  conditions  can 
be  supplied  by  the  near  field  analysis  of  Section  2  where  sub-generated 
turbulence  is  important.  Further  nimerical  evaluations  of  the  concentration 
distributions  for  typical  ocean  and  sub  conditions  will  be  presented  in 
the  final  report  of  this  study. 
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1 .  INTRODUCTION 

The  main  objective  of  this  study  is  to  develop  an  analytical  model 
for  the  prediction  of  the  mixing  and  decay  of  a  passive  scalar  within  the 
turbulent  wake  of  a  submarine  operating  beneath  the  ocean  surface.  Even 
though  the  interest  is  on  the  passive  scalar,  the  essence  of  the  problem 
is  the  determination  of  the  dynamics  of  a  submarine  wake.  The  decay  of  the 
concentration  of  a  passive  scalar  could  then  be  easily  obtained  from  the 
conservation  equation  for  a  known  flow  field. 

Numerous  theoretical  and  experimental  studies  on  the  wakes  behind 
unpropelled  bodies  in  a  uniform  medium  have  been  reported.  The  knowledge 
acquired  through  these  studies  is  not  readily  applicable  to  the  present 
case  of  interest  because  of  the  following  two  characteristics  of  a  submarine 
wake: 

1)  For  a  submarine  moving  at  a  constant  velocity,  there  is  no 
net  axial  momentum  being  deposited  in  the  wake.  The  phenomena 
associated  with  these  momentumless  wakes  are,  in  many  respects, 
different  from  a  wake  having  a  momentum  deficit. 

2)  The  ocean  environment  is  generally  stably  stratified.  Because 
of  the  turbulent  mixing  in  the  wake,  fluid  particles  within 

the  wake  are  out  of  equilibrium  with  the  stratified  environment. 
This  results  in  an  increase  ot  the  potential  energy  of  the 
wake.  Gravitational  restoring  forces  continuously  convert  part 
of  this  potential  energy  into  kinetic  energy  and  the  wake  tends 
to  collapse  vertically  and  spread  out  horizontally. 

These  characteristics  of  a  submarine  wake  have  not  been  subjected  to 
3ny  rigorous  theoretical  treatment  to  date.  In  a  recent  report,  Ko^  has 
presented  a  simplified  phenomenological  model  with  the  main  emphasis  on 
identifying  the  important  physical  mechanisms  which  affect  the  wake 
dynamics.  Even  though  a  fairly  good  agreement  with  some  existing  laboratory 
data  have  been  obtained,  this  simplified  model  must  be  improved  before 
it  can  be  used  as  a  reliable  prediction  tool.  This  simplified  model  along 
with  several  improvements  is  described  in  Section  2  of  this  report. 

Another  related  problem  having  to  do  with  the  effect  of  the  ocean 
background  turbulence  on  the  wake  diffusion  has  never  been  addressed 
previously.  Based  on  the  previous  calculations  (see  Ko^),  the  turbulence 
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generated  by  the  submarine  and  its  propulsion  system  persists  only  for  a 
relatively  short  distance  because  of  the  stabilizing  effect  of  the 
stratified  environment.  Therefore,  we  might  expect  the  far  wake  diffusion 
to  be  dominated  by  the  background  turbulence.  However,  our  present  level 
of  understanding  in  this  regard  is  rather  limited.  As  a  first  step,  we 
will  present  in  Section  3  a  partial  survey  of  the  existing  understanding 
of  the  oceanic  turbulent  diffusion.  The  application  to  the  submarine  wake 
will  also  be  indicated. 
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2.  SUBMARINE  WAKE  MODEL 

The  basic  approach  follows  closely  the  analysis  of  Reference  1  with 
a  few  improvements.  Some  of  the  assumptions  in  the  previous  analysis  will 
be  relaxed  in  order  to  achieve  a  more  reliable  first-order  predictive  tool, 
at  least  for  the  case  of  a  negligible  background  turbulent  diffusion. 
Therefore,  most  of  the  details  in  the  derivation  of  the  governing  equations 
will  be  omitted  here  and  reference  may  be  made  to  the  previous  analysis. 

For  a  submarine  moving  with  uniform  speed  UQ  at  a  depth  H  below 
the  ocean  surface,  a  Cartesian  sub-fixed  coordinate  system  is  chosen  as 
shown  in  Figure  1.  The  undisturbed  medium  outside  the  wake  is  assumed  to 
be  horizontally  stratified  with  the  density  given  by 

PQ  =  p(l  -  az)  (2-1) 

with  p”  denoting  the  mean  density  at  the  submarine  depth  and  a  being  the 
given  constant  gradient.  The  density  field  inside  the  wake  is  taken  to  be 

Pi  =  F(1  -  Bz)  (2-2) 

with  b,  the  internal  density  gradient,  reflecting  the  degree  of  mixing 
inside  the  wake.  The  effects  of  density  discontinuity  that  exists  at  the 
wake  boundary  for  this  assumed  density  distribution  will  be  considered  in 
a  later  study.  One  of  the  main  improvements  of  the  present  modeling  will 
be  the  inclusion  of  a  suitable  mixing  equation  to  account  for  the  variation 
of  B. 

The  concepts  of  local  similarity  and  turbulent  entrainment  of  the 
external  fluid  will  be  used  cs  before.  Assuming  that  there  is  no  direct 
effect  of  turbulence  on  the  return  of  a  displaced  fluid  pa 'tide  inside  the 
wake  to  its  equilibrium  position,  the  .'blowing  governing  equation  is 
obtained  from  the  momentum  equations 

uo  (1  +  TT  =  (“  •  5)9  -  (1  -  “I*  f>* 
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where  fj  is  the  local  velocity  scaling  factor  such  that 

v  =  yf2  and  w  =  -zf2  (2-4) 

The  values  of  a  and  b  represent  the  major  and  minor  axis  of  the  elliptical 
wake  boundary  given  by 


u2  7  2 
*—  +  =  1 
a’  b2 


(2-5) 


The  governing  equations  for  a  and  b  are  given  by 


uo37=af>  +  Ey 

“o  £r  -  -af>  *  E; 


(2-6) 


with  and  E2  denoting  the  turbulent  entrainment  velocities  in  the  y  and  z 
directions  at  the  wake  boundary  y  =  a  and  z  =  b.  If  the  entrainment 
velocity  is  assumed  to  be  proportional  to  the  distance  from  the  center  of 
the  wake,  then  we  have 


bF  b^JodA 
z  ~  a  y  2  A  dX 


(2-7) 


where  A  =  wab  denotes  the  area  of  the  wake.  The  entrainment  velocity  is 
then  related  to  the  local  turbulent  intensity  and  scale  using  the  experimental 
results  of  Naudascher,  which  gives 


U 


_c  dA 
A  dX 


Km 


m 


(2-8) 


with  Kx  w  2.8. 

The  formulation  presented  thus  far  is  identical  to  that  given  in 
Reference  1.  The  intent  of  what  follows  here  is  to  more  accurately  account 
for  the  effect  of  density  stratification  on  the  rate-of-change  of  the 
turbulent  energy  and  its  effect  on  the  wake  density  gradient  s. 
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The  total  integrated  turbulent  energy  in  the  wake  is  given  by 


E 


t 


=  p  u, 


•2 


m 


(2-9) 


where  is  the  averaged  turbulent  velocity  over  the  wake.  An  equation  for 
the  rate-of-change  of  the  turbulent  energy  within  the  wake  in  a  stratified 
medium  can  be  written  as 


dEt 

uo  TT  =  p  •  D  -  s 


(2-10) 


where  P  =  rate  of  production 
D  =  rate  of  dissipation 

S  =  stratification  effect  on  turbulent  energy  exchange 

The  modeling  of  the  production  and  the  dissipation  terms  follows  the  same 
line  of  reasoning  as  Reference  1  in  which  it  is  shown 

(2-11) 

(2-12) 

The  constants  are  determined  from  the  experiment  of  Naudascher  to  be  K2  =  140 
and  C  =  8.  The  modeling  for  the  stratification  effect  on  the  turbulent 
energy  balance  marks  the  main  difference  of  the  present  analysis  from  the 
previous  one.  The  previous  analysis  was  based  on  an  overall  potential  and 
kinetic  energy  conservation  concept  which  led  to  a  representation  of  the  term 

//  p  V  g  dA 

which  appears  in  the  turbulent  energy  equation.  However,  the  fact  that 
the  effect  of  stratification  on  turbulence  vanishes  when  8  =  o  leaves 
some  doubt  as  to  its  validity.  Furthermore,  because  of  the  need  to  model 
the  (p'w')  term  for  the  mixing  equation,  there  is  little  reason  not  to  mode! 
the  effect  of  stratification  on  the  turbulent  kinetic  energy  directly. 


Um 

P  -  K2  upa  — =- 
Uo2 

and 

D  =  Cirp  u;3  a 
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Since  this  modeling  was  not  considered  in  the  previous  report,  we 
will  go  into  some  detail  in  its  formulation  here. 

The  modeling  is  based  on  the  concepts  of  similarity  and  eddy 
diffusivity.  Let's  sssune 

z*>  (2-13) 


where  y*  =  y/a,  z*  =  z/b,  and  Fc  is  a  distribution  function  which  vanishes 
at  the  wake  boundary  and  has  a  maximum  at  the  center  of  the  wake.  One 
simple  form  of  Fc  can  bs  taken  as 


F 


c 


1  -  y*2  -  z*2 


Using  an  eddy  diffusivity  formulation,  we  let 


max 


=  e  6  p 

=  K'u'bSfT 
m 


(2-14) 


(2-15) 


The  value  of  K'  can  be  estimated  from  the  diffusivity  for  an  incompressible 
turbulent  wake  (e.g.,  Townsend)  to  be  about  0.2.  With  the  above  assumptions, 
evaluation  of  the  term  S  is 


S  =  9  /A/^  dA  =  4  g  K'u^bp  dz/^1  "  Z*2  (1  -  y*2  -  z*2)dy 


=  £  ab2 


(2-16) 


Therefore,  the  turbulent  energy  equation  is  then 


d(u'2ab)  u'5 

uo  -t— =  Kza  v  ■  Cu«3a  -  t 


(2-17) 
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The  proper  mixing  model  for  determining  the  local  density  gradient 
within  the  turbulent  wake  is  obtained  from  the  averaged  density 
conservation  equation,  which  gives 


3p.;  _  3p,  _  3p^ 

U  — +  v  — -  +  ■ 
o  ax 


a  — rrrr  ,  a  — rrrr  .  a  — 

u  +  v  +  ilp 


v] 


(2-18) 


Assuming  ~  p'u*  is  small  and  noting  that 

3X 


=  p(l  -  Bz),  v  =  yf;,  w  =  -zfx 


an  integration  of  the  LHS  of  equation  (2-18)  over  the  upper  half  of  the 
wake  gives 

-  vb  (2-19) 

where  a*  =  ab  and  B*  =  Bb.  The  integration  of  the  RHS  of  equation  (2-18) 
results  in  only  one  non-vanishing  term 

T  =  f3  p  V  (z  -  0)  dy  (2-20) 

which  represents  the  transport  of  mass  across  the  interface  z  =  0.  With 
equations  (2-13)  -  (2-15),  this  transport  term  can  be  readily  evaluated 
to  be 


T  =  f  K'u^  P  Bab 


(2-21) 


Then,  by  equating  (2-19)  and  (2-21),  we  obtain  the  equation  for  mixing 


U  E 

(°  -  «  >  ib  W  *  *  T  •  uo  ~5T  -  2K  V 


(2-22) 
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2.1  NON-DIMENSIONAL  GOVERNING  EQUATIONS 

By  using  the  following  characteristic  quantities 

Length  ...  L  (characteristic  dimension  of  the  flow) 
Velocity  ...  UQ  or  Ure<r  =  u^Q  (initial  turbulent  intensity) 

Time  ...  N"1  (N  =  /ag  ..  Vaisala-Brunt  frequency) 

and  let 


a*  =  a/L,  b*  =  b/L 

u  *  =  u'/U  x,  E  *  =  E  /U  x,  E  *  =  E  /U  x 
m  nr  ref’  2  z  ref*  y  y7  ref 

fj*  =  fj/N.  X*  =  XN/UQ 


and  define  two  non-dimensional  numbers 


F 


r 


NL 


...  tu’-bulent  Froude  number 


"n  -  WU0 


Then,  the  governing  equations  can  be  summarized  as 


a£  -  a*f>* +  FrEy* 
db* 


Ki  b* 
with  Ez*  =  T  u/  ^ 


Ki 

cy  2  um 


dx*  -  -b*fl*  +  FrV 


(1  +  eii  .  «1)  f,« 

U  a*2  dx*  U  a*J  a*2  U  a *2'  T1 


U  *5 

d  m  *2  ~  f  r  u  ‘i;,  — _ c 

dx*  um  hr[  N  N2  b*  L 


% 


*3 


-  K, 


u  *3 
m 


b*  1X1  a* 


•  F7  T  S  um*b*] 
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l^rT^  [!£*>“*-  <2K' +  K»  £>  **] 

The  fact  that  the  turbulent  energy  equation  allows  a  non-vanishing 

derivative  of  u  *  when  u *  =  0  requires  a  constraint  to  the  numerical 
mm 

program  such  that  u^*  is  set  to  be  identically  zero  once  it  reaches 
zero.  Another  somewhat  artificial  constraint  is  placed  on  the  magnitude 
of  S  such  that  it  will  not  exceed  a. 

Now  assuming  the  wake  to  be  initially  circular  and  taking  the 
characteristic  length  L  to  be  equal  to  the  initial  wake  radius,  the 
proper  set  of  the  initial  conditions  is 

ao*  =  bo*  =  1 
*io*  =  0 


and  a  specified  bq*. 

The  equations  point  out  that  the  proper  scaling  parameters  for 

ui. 

the  submarine  wake  are  the  turbulent  Froude  number  Fr  «  and  the 
initial  turbulent  intensity  level  In  addition,  the  initial 

degree  of  mixing,  represented  by  B0/a  appearing  in  the  form  of  an  initial 
condition,  constitutes  an  additional  parameter  in  the  problem.  The 
value  of  these  parameters  for  any  particular  experiment  or  field 
condition  is  not  apparent  at  the  present  time,  tven  though  these  parameters 
are  not  completely  independent  of  one  another,  the  actual  determination 
of  their  magnitudes  will  have  to  rely  quite  heavily  on  empirical  results. 
While  searching  for  the  proper  way  of  choosing  these  parameters, 
parametric  studies  will  be  performed  and  the  results  will  be  presented 
in  the  final  report. 
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3.  PRELIMINARY  CONSIDERATIONS  OF  TURBULENT  DIFFUSION 
OF  A  PASSIVE  SCALAR  IN  THE  OCEAN 


To  describe  the  diffusion  of  a  passive  scalar  (e.g.,  a  radioactive 
trace  contaminant)  in  the  ocean,  one  has  to  quantitatively  understand 
the  kinematics  of  ocean  turbulence.  If  one  seeks  to  find  the  mean 
concentration  of  a  contaminant,  F.  (x,y,z),  diffusing  from  a  fixed  source 
in  an  ocean  flow  field,  the  usual  approach  is  to  seek  solutions  of  the 
time  averaged  species  conservation  equation 


(3-1) 


where  c!ul  represents  the  species  flux  in  the  j  direction  and  «.  represents 

•  J  I 

the  time  averaged  rate  of  production  (or  loss)  of  species  i  by  chemical 
reaction  or  radioactive  decay.  ^  represents  a  material  derivative. 

The  derivation  of  Eq.  (3-1)  assumes  that  the  nature  of  the  flow 
field  is  described  by  a  stationary  random  function  of  time.  That  is, 
when  the  turbulence  properties  are  averaged  over  a  suitable  time  increment 
T,  a  stationary  value  for  the  turbulence  properties  is  obtained.  In  a 
meandering  ocean  flow  field,  various  averaged  turbulence  properties  may 
result  depending  on  the  magnitude  of  the  time  tcale  T.  For  example,  large 
scale  horizontal  ocean  diffusion  properties  only  achieve  a  stationary  value 
when  averages  are  taken  over  periods  of  several  hours  to  a  day. 

For  the  analysis  to  be  presented  in  this  report,  stationary  values 
of  the  turbulence  properties  will  be  assumed  in  order  to  provide  a  basic 
estimate  of  the  diffusior.al  characteristics  of  oceanic  turbulence. 


3.1  DIFFUSION  COEFFICIENTS 

To  solve  the  species  conservation  equation,  Eu.  (3-1),  expressions 
for  the  species  flux  cTuT  must  be  adopted.  Using  the  Boussinesq  approximation, 

*  J 

turbulent  diffusion  coefficients  can  be  defined  by  the  equation 

Kj  =  -  c!uj/(ac./3Xj) 
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These  coefficients  do  not  necessarily  assume  a  constant  value  and  may  be 
functions  of  the  distance  x  from  a  source  (or  time  t  =  x/U).  Specific 
forms  K  =  K(x)  have  been  suggested  by  Sutton*  '  for  various  stability 
conditions  corresponding  to  the  atmospheric  diffusion  of  smoke  plumes. 

The  change  of  the  diffusion  constant  with  time  is  best  calculated  by 
Taylor's^  diffusion  theory  of  continuous  movements  in  a  homogeneous 
turbulent  field. 

Taylor  defined  the  function  R(t)  as  the  coefficient  of  correlation 
between  the  turbulent  velocity  ut  of  a  particle  of  fluid  at  time  t  and 
the  velocity  u^+t  of  the  name  particle  at  a  time  t  later,  such  that, 

r(t)  =  utut+x/u2  (3-3) 

where  u2  represents  the  mean  square  velocity  fluctuation  intensity  at  time 
t.  If  5  is  the  displacement  of  a  particle  of  fluid  at  time  t  from  its 
original  position  at  time  t  =  0,  Taylor  proved  that 

-t  rf 

C7 :  a2  =  2  u7  J  /  R(x)dxdt'  (3-4) 

o  o 

For  small  values  of  t,  such  that  R(t)  =  1.0 


$2  =  a2  =  u2t2  (3-5) 

For  times  t  »  t*,  (where  t*  denotes  the  Lagrangian  time  scale  =  J  R(x)dx) 

*o 

the  correlation  R(x)  is  effectively  zero.  Hence  if  diffusion  proceeds  for 
times  much  greater  than  t*,  the  integral  in  Eq.  (3-4)  will  reach  a 
limiting  value  I.  In  that  case,  the  mean  square  displacement  is  proportional 
to  the  time  of  diffusion,  i.e.. 


o2  =  2  i'2  It 

or  rT=r  (3-6) 

a2  =  2  V  u2  A,  t 


where  A.  =  V  u2  / 

L  o 

If  the  particles  had 

diffusion  equation 


R(x)dx  is  defined  as  the  lagrangian  integral  length, 
been  dispersed  by  a  diffusion  process  obeying  the 
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5c  _  k  52c 


the  solution  of  Eq.  (3-7)  for  the  concentration  c  would  be 


(3-7) 


c  =  const  exp  (_y2/ 4K  t)  (3-8a) 

/t  y 

or 

(3-Sb) 

Equations  (3-8a)  and  (3-8b)  yield  a  mean  square  displacement  ay  of  fluid 
particles  in  the  y  dimension;  given  by 


oy2  =  2  Kyt  (3-9) 

Based  on  this  relation,  one  can  define  a  diffusion  coefficient  from 
Eq.  (3-9)  in  terms  of  the  mean  square  displacement,  i.e.. 


The  diffusion  coefficient  uy  given  by  Eq.  (3-10)  can  thus  be  used  directly 
in  a  simple  diffusion  equation  of  the  form  given  by  Eq.  (3-7). 

(4l 

Batchelor'  '  extended  Taylor's  ideas  with  the  aid  of  dimensional 
arguments.  Ke  obtained  expressions  for  the  diffusion  coefficients  K 
(or  da2/dt)  applicable  to  the  horizontal  diffusion  of  a  passive  scalar 
(or  trace  contaminant)  in  a  homogeneous  ocean.  The  characteristic  velocity 
scale  is  related  to  the  rate  of  turbulent  energy  dissipation  per  unit 
mass  e,  and  to  the  standard  deviation  of  the  initial  source  field  aQ  by 
the  expression 


(3-11) 


Batchelor  finds  three  regimes  of  relative  diffusion,  the  two  given  by 
Taylor's  relations  [Eq.  (3-5)  and  (3-6)]  and  one  important  intermediate 
regime  given  by  the  following  rates  of  growth 
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K  =  I^W=  Ci^q)^3  (initial)  t  «  t*  (3-12a) 

K  =  ^  =  C2e(t  -  tj)2  (intermediate)  t  ^  t*  (3-1 2b) 

with  t2  =  C30o2/3£'i/3 

K  =  C4^L  (asymptotic)  t  »  t*  (3-1 2c) 

For  the  intermediate  times  as  given  by  Eq.  (3-12!>) 

K  <V  et2  -V,  e1/3a4/3  (3-13) 

and 

o2  *v»  et3  (3-14) 


Thus  the  diffusivity  grows  as  the  4/3  power  of  the  particle  field  (or 
plume)  size,  and  the  plume  grows  as  the  3/2  power  of  the  time.  The  4/3  law 
was  first  derived  by  Richardson^  in  1926.  Richardson’s  law  of  relative 
diffusion  indicates  that  diffusion  is  an  accelerating  process,  the  rate 
of  growth  increasing  with  the  size  of  the  field,  at  least  while  the  supply 
of  eddies  of  the  requisite  size  lasts. 

3.2  TWO  DIMENSIONAL  HORIZONTAL  DIFFUSION 

Richardson's  law  has  been  found  to  describe  fairly  well  the 
horizontal  diffusion  in  the  ocean.  Experimentally,  horizontal  diffusion 
has  been  studied  extensively  using  the  fluorescent  dye  technique.  An 
overall  discussion  of  such  experiments  is  presented  in  the  review  article 
of  Bowden^.  In  some  cases  the  dye  has  been  released  continuously  from 
a  point  source  into  a  current  and  the  dispersion  of  the  plume  observed. 

In  other  cases  the  dye  has  been  released  as  instantaneously  as  possible 
and  the  dispersion  of  the  resulting  patch  observed  as  a  function  of 
horizontal  coordinates  and  time.  In  both  cases  the  initial  dispersion 
is  usually  three  dimensional  but  after  a  relatively  short  time  the  vertical 
dispersion  becomes  severely  restricted  by  the  presence  of  a  thermocline, 
and  the  subsequent  spreading  is  effectively  two  dimensional.  In  that 
case,  the  diffusion  equation  for  a  point  source  in  a  uniform  turbulent 
flow  of  velocity  U,  assumes  the  form 
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U 


(3-15) 


For  the  case  where  an  isolated  source  emits  a  mass  per  unit  time  of  species  i 
at  a  constant  rate  Q.*  per  unit  distance  in  the  vertical  direction, 

Eq,  (3-15)  has  the  solution 


ci 


QjVp 

*'/2iT(TyU 


(3-16) 


where 


1  V 

2  “3t  ■ 


(3-17) 


( 7) 

0kubov  '  in  1968  has  made  a  review  of  20  sets  of  data  on  instantaneous 

releases  of  dye  in  the  ocean  over  the  last  ten  years.  They  include 

experiments  in  the  North  Sea  and  in  the  Cape  Kennedy  area,  and  off  of 

Southern  California.  Okubo  has  tabulated  the  variance  a 2.  and  the 

rc 

apparent  diffusivity  K  as  a  function  of  time  t  for  each  of  the  20  releases. 

a 

These  results  are  shown  in  Figures  2  and  3  where  is  shown  plotted 
against  the  diffusion  time  and  K,  is  plotted  against  the  scale  of  diffusion, 

a 

i  =  3ar(;.  These  diagrams  cover  time  scales  ranging  from  1  hr.  to  1  month 
and  length  scales  from  100m  to  100km.  For  the  increase  of  variance  with 
time,  Okubo  found  that  a  best  fit  of  the  data  could  be  obtained  if 


arc  =  -0108t2’34  (3-18) 

This  is  a  somewhat  smaller  rate  of  growth  than  given  by  the  t3  law  based  on 
the  Richardson  theory,  Eq.  (3-14).  The  relation  between  the  apparent 
diffusivity  and  scale  was  found  by  Okubo  to  be 

K  =  .0U1*15  ( 3-1 9a) 

Q 

4/3 

compared  with  the  theoretical  1  law.  However,  the  differences  with 
the  Richardson  theory  are  not  large. 


paHikMHfeiliiaijHi 
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If  one  assumes  that  the  horizontal  diffusivity  can  be  written  in 
the  form 


and  that 


K  =  at1 

a 


£  =  ba 


[  3-1 9b) 


then  by  Eq.  (3-10)  one  can  readily  show  that 


=  [(2  -  n)St]1/(2-n) 


(3-20) 


where 


^  Ln, 
a  =  b  a 


The  concentration  distribution  of  a  passive  scalar  downstream  of  a  continuous 
point  source,  moving  at  a  velocity  U,  in  a  quiescent  ocean,  will  be  given 
by  Eqs.  (3-16)  and  (3-20),  i.e.. 


Q.*/pU 


'i  fa  [(2  -  n)ax/u] 


i/(2-n) 


(2  -  n)ax/ul2/  2-n^ 


For  the  constants  n  =  1.15  the  concentration  drops  off  rapidly  with 
distance  from  the  source  x  as 


(3-21) 


cT*-i- 
1  x1'17 


(3-22) 


and  the  size  of  the  field  increases  as 


Oy  *  X 


(3-23) 


If  the  Richardson's  4/3  law  is  assumed  to  be  valid,  n  =  4/3,  hence 


ci  "Jn  ’  °y'  x‘'5 


(3-24) 


Calculations  of  concentration  versus  distance  behind  the  source  are 
presented  in  Fig.  4  for  both  the  n  =  1.15  and  n  =  4/3  laws. 
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In  tne  following  section  more  general  solutions  are  presented  of  the 

diffusion  equation  including  the  effects  of  vertical  diffusion,  and 

radioactive  decay.  In  addition  solutions  will  be  presented  which  allow 

for  arbitrary  initial  conditions  based  on  the  near  field  submarine  wake 
'll 

solutions  of  Kg'  in  a  stratified  medium.  Diffusion  calculations  for 
typical  submarine  wakes  will  be  included  in  the  final  report  of  this 
study. 


3.3  THREE  DIMENSIONAL  DIFFUSION  WIiH  DECAY 

The  semi -empirical  steady  state  diffusion  equation  is  written  below 
for  a  three  dimensional  Cartesian  coordinate  system  x,  y,  z  wit!,  stream 
velocity  U  in  the  x  direction 


3c-  .  3c-  ’  3C-  c. 

u  V (Ky  +  «  (Kz lf>  - ^ 


(3-25) 


Here  the  z  axis  is  directed  vertically  downward,  Ky  and  Kz  are  the  turbulent 
diffusivities  defined  by  Eq.  (3-2)  in  the  y  and  z  directions.*  To  take 
into  account  a  possible  exponential  decrease  of  cT  with  characteristic 
decay  time  tc  (e.g.,  caused  by  radioactive  decay),  the  last  term  on  the 
right  of  Eq.  (3-25)  has  been  chosen  to  represent  the  production  term, 
in  Eq.  (3-1). 

A  special  solution  of  the  diffusion  equation  (3-25)  is  possible  for 
an  isolated  fixed  source  of  strength  (mass  per  unit  time)  located  at 
(0,  0,  h)  [a  distance  h  below  the  water  surface  as  shown  below] 


*A  term  of  the  form  ~  (Kx  -^-)  is  neglected  in  Eq.  (3-25)  in  comparison  to 

3C  - 

U  .  This  assumption  is  quite  good  for  problems  involving  continuous 
ax 

sources . 
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For  the  case  of  a  uniform  stream  velocity  U,  and  for  values  of 
and  Kz,  which  are  constants  or  functions  of  x,  the  solution  of  Eq.  (3-25) 
is  found  to  be^ 


(x-  2)  =  “p  [•  i^]exp  [-  i  (t>2] 

y  z  j 

.|exp[-l(^)2].exp  [-i(^)2]j  (3-26) 


The  dispersion  coefficients  are  <jy  =  (— jjr-)  and  o2  =  (— jf~)  (3-27) 

This  solution  corresponds  to  the  following  initial  and  boundary  conditions 
imposed  on  Eq.  (3-25). 


x  =  0: 

c-  -*•  Q. 

i 

z  =  0: 

1) 

O 

z  -*■  «: 

c-  -*■  0 

If  the  source  is  far  below  the  surface  h  »  az,  then  Eq.  (3-26)  takes  the 
simple  Gaussian  form 

-*<£*]  <«8> 


where  z  =  z  -  h. 

3.4  VERTICAL  DIFFUSION 

In  order  to  make  quantitative  use  of  Eq.  (3-26),  in  the  calculation 
of  three  dimensional  concentration  distributions,  the  diffusion  coefficients 
Ky,  and  Kz  must  be  given  or  the  dependence  of  and  oz  on  the  distance  x 
must  be  specified.  In  the  previous  section  it  was  shown  that  the  horizontal 
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diffusion  coefficient  Ky  in  the  ocean  could  be  characterized  by  the 
Richardson  4/3  law  or  a  closely  related  empirical  expression  of  the  form 
Ky  =  atn.  Typical  ocean  measurements  (see  Fig.  3)  yield  values  of  Ky  in 
the  range 


Ky  =  0(103  to  107)  f^- 

Corresponding  values  of  the  vertical  diffusion  coefficient  Kz  are  usually 
much  lower  and  are  reported  to  lie  in  the  range^ 

Kz  =  0(1  to  103) 

The  principal  mechanism  by  which  turbulence  is  generated  in  the  upper 
layers  of  the  ocean  is  by  the  action  of  the  wind  shear  acting  on  the 
surface  layer.  From  observational  data,  Sverdrup^ 10^  proposed  the 
following  empirical  relationship  between  NZq  [the  surface  vertical 
turbulent  eddy  stress  coefficient  =  -  u'w'/fau/az)]  and  the  wind  speed 
w(m/sec); 


pN*0  =  *.3w2  f°r  *  >  6  m/sec  (3-29) 

For  most  studies,  Reynolds'  analogy  between  momentum  and  diffusive 
transport  is  assumed,  i.e.. 


kzq  nzo 

In  the  case  of  a  stably  stratified  ocean  (the  density  increasing  with 
depth)  the  mean  buoyancy  field  acts  to  suppress  the  generation  of  turbulence. 
The  criteria  used  to  determine  the  fluid  conditions  for  which  turbulence 
is  completely  suppressed  is  based  on  the  Richardson  number 


r  =  9  8p/3z 

1  P  (3U/3Z)2 


(3-30) 


where  p  is  the  density  of  the  fluid,  u  is  the  current  velocity  and  g  the 
acceleration  of  gravity  in  the  vertical  z  direction.  The  Richardson 
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number  represents  the  ratio  of  the  rate  of  suppression  of  turbulence  by 

the  buoyancy  field  to  the  rate  of  production  by  the  action  of  Reynolds 

stresses.  If  R*  is  greater  than  some  critical  Richardson  number,  R-j  .., 
1  fill  cnt 

(Ricrit  order  *’5)  then  the  theory  of  Ellison'  '  indicates  that  the 

diffusivity  Kz  should  vanish. 

f  12I 

Sundaram  and  Rehm'  *  in  their  studies  of  lake  thermoclines  have 
recently  adopted  the  following  empirical  expression  to  determine  the 
effect  of  stability  on  the  vertical  diffusivity  of  heat 


VKh0  =  (1  +  dR,)*1 


(3-31) 


where  ax  is  a  constant  found  to  be  approximately  .6.  In  Sundaram  ar,d 
Rehm's  analysis,  the  velocity  gradient  3u/3z  is  assumed  to  be  determined 
by  the  expression  for  a  constant  stress  turbulent  layer 


3u/»z  =  u*/kz  (3-32) 

where  k  is  Karman's  constant  -  .4  and  u*  =  »t^7p,  is  the  friction  velocity 
with  t  being  the  surface  shear  stress.  Equations  (3-31)  and  (3-32)  were 

^  (l3) 

shown  to  adequately  represent  both  the  field  measurements  in  Cayuga  Lake' 
and  the  laboratory  measu.  °ments  in  a  stratified  channel  flow  (see  Fig.  4). 

An  overall  empirical  expression  for  tht  vertical  diffusivity  coefficient 
Kz  can  be  obtained  by  combining  Eqs.  (3-29)  *  (3-32),  i.e., 

K  =  j= - - -i  (3-33) 

with  w  in  m/sec.  For  typical  values  of  the  stratified  density  gradient 
—  !§•(*  10"7g  cm"4)  and  shear  velocity  u*  (*  33  cm/sec),  encountered  in 
the  ocean,  it  is  found  from  Eq.  (3-33)  that  the  vertical  diffusivity  Kz 
is  reduced  by  a  factor  of  one-half  by  a  depth  z  wlOO  meters  (near  the 
top  of  the  ocean  thermocline).  At  a  depth  of  300  meters,  the  magnitude 
of  K  is  only  one  tenth  of  its  value  near  the  surface,  thus  demonstrating 

^  -5, 

that  a  density  gradient  as  small  as  10  g  cm  can  virtually  eliminate 
vertical  transport  in  the  thermocline. 
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( 14) 

Ozmidov'  reports  a  striking  observation  of  this  effect  in  his 
experiments  of  the  diffusion  of  patches  of  dye  (rhodamine  C)  in  the  Black 
Sea  during  the  fall  of  1964: 

“In  its  initial  phase,  the  diffusion  of  the  dye  from 
an  instantaneous  point  source  is  symnetric  in  all 
directions.  But  after  the  patch  reiches  dimensions 
of  several  meters  it  begins  to  "flatten  up”  sharply, 
and  whereas  tne  increase  of  the  patch  along  the 
vertical  has  practically  stopped,  the  patch  continues 
to  increase  rapidly  along  the  horizontal  directions. 

After  several  hours  had  elapsed,  the  patch  reached 
horizontal  dimensions  of  the  order  of  1  km,  whereas 
its  vertical  extent  did  not  usually  exceed  15  to  20  m 
(in  this  case  the  layer  of  discontinuity  was  located 
at  a  depth  of  M>0m).“ 

(l4) 

Ozmidov'  '  has  proposed  a  simple  theory  for  turbulent  diffusion 
in  the  ocean,  which  seeks  to  quantitatively  model  the  anisotropic  effects 
noted  ir.  his  ocean  dye  diffusion  experiment.  His  arguments  are  as  follows; 

If  the  energy  supply  of  tne  ocean  is  sufficiently  large,  an  interval 
of  scales  f-nn  to  £cr  will  exist  in  which  the  turbulence  is  isotropic. 
The  lower  bound  iQ  is  equivalent  to  the  smallest  dissipation  length  scale 
of  the  inertial  subrange,  i.e.. 


l0  =  (v3/e)1/‘*  (3-34) 

where  v  is  the  kinematic  viscosity  of  sea  water  and  e  is  the  rate  of 
turbulent  energy  dissipation.  The  upper  bound  for  which  the  turbulence 
will  be  isotropic  is  characterized  by  the  eddy  size  &cr  for  which  the 
stratification  density  barrier  is  just  overcome  Ucr  is  related  to  the 
Richardson  number  based  on  e),  i.e.. 


t 


cr 


(3-35) 


In  the  interval  tQ  <  t  <  the  diffusion  of  impurities  from  a  instantaneous 
point  source  is  identical  for  all  directions  and  a  cloud  of  diffused  mattt»* 
will  have  a  spherical  shape.  Eddies  with  dimensions  greater  than  tcr  can 
no  longer  take  part  in  the  vertical  diffusion  transfer,  thus  limiting  the 
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effective  vertical  diffusion  coefficient  Kz.  The  limiting  value  of  Kz 
can  be  estimated  from  the  4/3  power  law  for  isotropic  turbulence 


k  -  r  o1/3»  4/3  - 
Kzmax  "  Cie  *cr 


(3-36) 


where  Ci  is  a  constant  approximately  equal  to  .1. 

Ozmidov  argues  that  in  the  range  of  scales  from  i  to  a  in  which 

o  cr 

the  turbulence  is  isotropic  the  values  of  the  vertical  diffusion  coefficient 
K  and  of  the  horizontal  coefficient  K  are  equal  (and  are  given  by  the 

4  X> 

4/3  law).  However,  for  phenomena  with  scales  greater  than  a  ,  K  and 

vl  A> 

Kz  become  markedly  different  as  illustrated  in  the  sketch  below. 


laminar  v 


For  scales  a  >  acr,  Kz  reaches  its  maximum  value  Kzmax  (given  by  Eq.  (3-36)) 
determined  by  ouoyancy.  The  horizontal  coefficient,  K^,  jumps  to  a  new  two 
dimensional  4/3  law  curve  (with  constant  C2)  and  continues  to  increase 
with  a  up  to  some  maximum  scale  tj  which  may  be  of  the  order  of  several 
kilometers.  Sample  calculations  for  typical  values  of  3p/3z  and  e  in  the 


ocean,  yield  limiting  values  of  K-.  „  in  the  range  of  10  to  100  cm2/sec  and 

max  (14) 

values  of  acr  in  the  range  of  1  to  10  meters'  .  These  numbers  confirm 
the  relatively  small  extent  of  vertical  diffusion  in  a  stratified  ocean. 


One  essential  feature  pointed  out  by  this  model  is  that  while  far 
field  diffusion  effects  may  well  be  given  by  the  4/3  law  [e.g.,  Eq.  (3-24)] 
for  horizontal  diffusion  and  the  corrected  Sverdrup  law  [Eq.  (3-32)]  for 
vertical  diffusion,  these  models  are  not  adequate  in  the  near  field  where 
diffusion  changes  from  isotropic  to  anisotropic  behavior. 
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To  estimate  the  effect  of  the  near  field  behavior  on  the  far  field 
solution,  it  is  possible  to  seek  alternate  solutions  of  the  two  dimensional 
diffusion  equation  [Eq.  (3-15)],  assuming  that  the  initial  field  concentration 
and  size  are  prescribed  as  distributed  initial  conditions  (rather  than 
point  source  initial  conditions).  Appropriate  near  field  solutions  for  a 
collapsing  submarine  wake  have  recently  been  set  forth  by  Ko^  and 
modifications  presented  in  Section  2.  These  solutions  would  provide  the 
needed  initial  wake  dimensions  and  initial  concentration  levels. 

The  reformulated  two  dimensional  (vertical  transport  neglected) 
diffusion  equation  and  boundary  conditions  can  be  written  as 


,,  ac  _  a2c  c 
U  «  -  Ky  -  Tc 


(3-37a) 


with  boundary  conditions 

9.  t 

r  =  cn  at  x  =  0  -  y  <  v  <  y  (3-37b) 

(virtual  origin)  c 

c-»0asx+»,  c  +  0  3sy  +  ±®  for  all  x 

where  iQ  represents  the  initial  width  of  the  di fusion  layer  at  the  joining 
point  with  the  near  field  solution. 

The  solution  of  Eq.  (3-37)  will  follow  the  previous  analysis  of 
Brooks^15^  for  diffusion  of  sewage  effluents  in  an  ocean  current  and  the 
thermal  discharge  analysis  of  Sur.daram  et.  al.^16\ 

Equation  (3-37)  can  be  converted  to  a  simple  diffusion  equation 
without  a  decay  term  by  the  transformation 

c  =  c  exp  [-  x/(Uxc)]  (3-38) 

Then  Eq.  (3-37)  reduces  to  the  form 


U  |£-  = 

3X 


Vx)  $ 


(3-39) 
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It  will  be  assumed  as  in  Eq.  (3-1 9b)  that  the  horizontal  diffusivity  is 
related  to  the  characteristic  length  i(x)  by  the  relation 

Ky  =  atn(x) 

and  (3-40) 

t(x)  =  bo(x) 

where  a  is  the  standard  deviation  of  the  concentration  distribution 
defined  by 


a(x)2 


y^dy 


(3-41) 


I*  one  requires  that  t  =  iQ  at  x  =  0,  where  c  is  taken  to  be  uniform  and 
equals  to  cQ  in  the  domain  -«.q/2  <  y  <  *  / 2,  then  b  ■■  Z/J  *  3.4. 

The  solution  of  Eq.  (3-39)  couV.  be  written  down  imnediately  if 
were  a  constant,  independent  of  x.  However,  for  Ky  being  a  function  of  x, 
a  simple  transformation  of  the  x  coordinate  can  be  used  such  that 


K0dx '  =  Ky(x)dx 


=  a 


(3-42) 


Then  Eq.  (3-39)  reduces  to  the  well  known  diffusion  equation  with  constant 
coefficient  Kq  (=  a*0n),  i.e.. 


U 


3C  _ 


3X ' 


=  K 


32C 


0  3, 


(3-43) 


For  the  initial  conditions  given  by  Eq.  (3-37b),  the  closed  form  solution 
of  Eq.  (3-43)  given  by  Crank^17\  is 


erf 


2y  +  V  _  2y^V  ( 

4^T7u  4/kox'/u  ( 


(3-44) 
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The  integral  in  Eq.  (3-42)  which  is  used  to  relate  the  physical 
x  and  transformed  x'  variables,  cai.  now  be  evaluated  by  employing  the 
error  function  solution  given  in  Eq.  (3-44),  i.e.. 


rr*  rrr/-.  y2Jdy  »  1  +  2  -2 


K  x' 


£  *  >  £  3 

0  cO  0 


V  U 


(3-45) 


Hence  the  transformed  axial  dimension  x'  is  found  in  terms  of  x  by 
integrating  Eqs.  (3-42)  and  (3-45), 


x'  =  J 


9 *U  ir  _  ?b2K  -i2/(2"n) 


(3-46) 


For  large  values  of  x,  far  from  the  source,  and  for  n  =  4/3  (Richardson  law) 


x'  ^  x3  x/£„  »  1 


'inder  tun  r.cme  circumstances,  the  concentration  s-»l,-tior. .  Eq.  (3-^4). 
takes  the  limiting  form, 

c  n,  (x')-l/2  x/£q  »  1 

Hence,  the  point  source  behavior,  given  by  Eqs.  (3-21)  and  (3-24) 

"c  (x)-3/2 

is  recovered  from  the  solution,  Eq.  (3-44),  for  diffusion  from  a  finite 
near  wake  source. 


In  summary,  the  complete  solution  for  diffusion  from  a  finite 

of  dimension  and  concentrati 
o 

in  the  following  non-dimensional  form 


source  of  dimension  £Q  and  concentration  cQ,  with  decay  is  written 


t^exp[-^D]|erfl#'erff#i  (3-47) 
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x 


■ 


where  x  =  x/£Q, 


=  2Ll=  1 
4o  2b\ 

y  *  y/t0,  L0 


it 


l  +  (1 


-  f)  2b2K 


.»] 


2/2-n 


=  LW 


a£ 


n-i 


and 


v  =  = 

o  M 

o 


U 


In  the  final  report  of  this  investigation,  a  parametric  study  will  be 
performed  to  determine  the  influence  of  the  various  prominent  parameters 
in  Eq.  (3-47)  on  the  magnitude  and  variation  of  c.-'taminant  concentration 
behind  a  typical  submarine  in  an  ocean  environment. 
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Fig.  1.  Coordinates  and  Density  Field 
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Fig.  2.  Variance,  of c,  of  Dye  Distribution  against  Diffusion 

Time  (Data  from  1961  Onwards)  (after  Okubo,  1968,  Ref.  7) 
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n  =  4/3 


n  =  1. 15 


Diffusion  Coefficient 
K  =  at" 

Source  Strength 

Q*/pU  =  1.0  ft 


DIFFUSION  TIME 


Fig.  4c.  Centerline  Concentration  Behind  a  Point  Source  Moving  at  a 
Velocity  U  in  an  Ocean  Environment,  Eq.  (3-21)  (v  =  0) 


